!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
program data_manage
!==================================================================================================================================|

   use mpi
   implicit none

   include 'data_manage_precision.include'
   include 'data_manage_parameter.include'

!----------------------------------------------------------------------------------------------------------------------------------|
!  ro, gp, mp, te, beta, az, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgv, fggx, fggy, fggv, 
!  fmgx, fmgy, fmgv, fmtx, fmty, fmtz, fmtv
!----------------------------------------------------------------------------------------------------------------------------------|
   integer(4), parameter :: n_var_total = 31
   integer(4), parameter :: var_total_dim(n_var_total) = (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 3,         &
                                                            1, 1, 3, 1, 1, 3, 1, 1, 1, 3 /)
   character(128), parameter :: var_total_name_list(n_var_total) =                                                                 &
      (/ 'ro', 'gp', 'mp', 'te', 'beta', 'az', 'vx', 'vy', 'vz', 'vv', 'bx', 'by', 'bz', 'bv', 'jx', 'jy', 'jz', 'jv',             &
      'fgx', 'fgy', 'fgv', 'fggx', 'fggy', 'fggv', 'fmgx', 'fmgy', 'fmgv', 'fmtx', 'fmty', 'fmtz', 'fmtv' /)

   character(4) :: ti_char, cpui_char
   character(128) :: datavis_path, var_name_list, directory, datavis_name
   character(256) :: header
   character(128), allocatable :: var_name(:)

   integer(4) :: output_method, output_level, n_var, len_str, format, i_var, j_var, i_block
   integer(4) :: ti, ci, i, j, i_e, j_e, ip1, jp1, cpui, nx_tmp, ny_tmp, i_c, i_p, j_c, j_p
   integer(4) :: nx_no_bnd, nx_block, nx_total, ny_no_bnd, ny_block, ny_total, nx_lv, ny_lv
   integer(4) :: id, cpu_number, ng, ratio, n_levels, nc, nc_lv
   integer(4) :: nt, n_blocks, lv, n_dims, sys, ini_block_x, ini_block_y, div_clean_flag
   integer(4) :: datamap_unit, datavis_unit, datavtu_unit
   integer(4) :: error, cpu_id, cpu_n, n_out
   integer(4) :: output_time_start, output_time_end
   integer(4), dimension(:), allocatable :: p_pos, n_blocks_all
   integer(4), allocatable :: var_dim(:), level_len(:)
   integer(4) :: data_offset

   logical(4) :: var_flag(n_var_total)

   real(8) :: gm, rangex, rangey, zero_posx, zero_posy, cfl, gmm1
   real(8), dimension(:), allocatable :: t, dt
   real(output_real8_flag), dimension(:, :), allocatable :: ro, mx, my, mz, bx, by, bz, en, gp, az, gx, gy
   real(output_real8_flag), dimension(:, :), allocatable :: roh, gph, mph, teh, betah, azh, vxh, vyh, vzh, bxh, byh, bzh,          &
                                                            jxh, jyh, jzh, fgxh, fgyh, fggxh, fggyh,                               &
                                                            fmgxh, fmgyh, fmtxh, fmtyh, fmtzh
   real(output_real8_flag), dimension(:, :), allocatable :: ro_lv, gp_lv, mp_lv, te_lv, beta_lv, az_lv,                            &
                                                            vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv,                              &
                                                            jx_lv, jy_lv, jz_lv, fgx_lv, fgy_lv,                                   &
                                                            fggx_lv, fggy_lv, fmgx_lv, fmgy_lv,                                    &
                                                            fmtx_lv, fmty_lv, fmtz_lv

   real(output_real8_flag), allocatable :: var(:, :, :, :), var_max(:), var_min(:), var_max_g(:), var_min_g(:)

   real(8) :: dx, dy, eps, dx_lv, dy_lv
   real(8), dimension(:, :), allocatable :: x_arr, y_arr
   real(8), dimension(:), allocatable :: x, y, xh, yh, x_lv, y_lv

   namelist /datavis_options/                                                                                                      &
      output_method,                                                                                                               &
      output_level,                                                                                                                &
      output_time_start,                                                                                                           &
      output_time_end,                                                                                                             &
      datavis_path,                                                                                                                &
      var_name_list,                                                                                                               &
      format,                                                                                                                      &
      header

!----------------------------------------------------------------------------------------------------------------------------------|
!  start
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_init(error)
   call mpi_comm_size(mpi_comm_world, cpu_n, error)
   call mpi_comm_rank(mpi_comm_world, cpu_id, error)

   if (cpu_id .eq. 0) then
      write(*, *) 
      write(*, '(a)') '|------------------|'
      write(*, '(a)') ' Data manage starts:'
      write(*, '(a)') '|------------------|'
      write(*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  check no reduplicate variable name
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, n_var_total
      i_var = 0
      do i = 1, n_var_total
         if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0) i_var = i_var + 1
      enddo
      if (i_var .gt. 1) then
         write(*, '(a, 1x, a)') ' ERROR: reduplicate variable name:', var_total_name_list(j)
         do i = 1, n_var_total
            if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0)                                          &
               write(*, '(1x, a, 1x, a, 1x, a)') trim(var_total_name_list(j)), '<-', trim(var_total_name_list(i))
         enddo
         write(*, *)
         stop
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  prepare
!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12
   datamap_unit = 10
   datavis_unit = 11
   datavtu_unit = 12

   open(datamap_unit, file = trim(parameters_filename))
   read(datamap_unit, nml = datavis_options)
   close(datamap_unit)

   call system("mkdir -p " // trim(datavis_path))

   var_flag = .false.

   len_str = len(var_name_list)
   do i = 1, len_str
      if (var_name_list(i:i) .ge. 'A' .and. var_name_list(i:i) .le. 'Z')                                                           &
         var_name_list(i:i) = char(ichar(var_name_list(i:i)) - ichar('A') +  ichar('a'))
   enddo

   n_var = 0
   do i_var = 1, n_var_total
      if (index(var_name_list, trim(var_total_name_list(i_var))) .ne. 0) then
         var_flag(i_var) = .true.
         n_var = n_var + 1
      endif
   enddo

   if (cpu_id .eq. 0) then
      write(*, '(a, i4)') ' Total number of variables is:', n_var
      write(*, '(a)') ' The following variables will be outputted for visualization:'
      write(*, *)
   endif

   allocate (var_name(n_var), var_dim(n_var), var_max(n_var), var_min(n_var), var_max_g(n_var), var_min_g(n_var))

   if (cpu_id .eq. 0) then
      do i_var = 1, n_var_total
         if (var_flag(i_var)) write(*, '(1x, a8, 1x)', advance = 'no') trim(var_total_name_list(i_var))
      enddo
      write(*, *)
   endif

   var_name = ''
   var_dim = 0

   j_var = 1
   do i_var = 1, n_var_total
      if (var_flag(i_var)) then
         var_name(j_var) = var_total_name_list(i_var)
         var_dim(j_var) = var_total_dim(i_var)
         j_var = j_var + 1
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  readparameters
!----------------------------------------------------------------------------------------------------------------------------------|
   open(datamap_unit, file = trim(datamap_path) // 'params.mapout')
   read(datamap_unit, '(a)') directory
   read(datamap_unit, *) cpu_number
   read(datamap_unit, *) n_dims
   read(datamap_unit, *) sys
   read(datamap_unit, *) ng
   read(datamap_unit, *) ini_block_x
   read(datamap_unit, *) ini_block_y
   read(datamap_unit, *) nx_no_bnd
   read(datamap_unit, *) nx_block
   read(datamap_unit, *) nx_total
   read(datamap_unit, *) ny_no_bnd
   read(datamap_unit, *) ny_block
   read(datamap_unit, *) ny_total
   read(datamap_unit, *) ratio
   read(datamap_unit, *) n_levels
   read(datamap_unit, *) div_clean_flag
   read(datamap_unit, *) gm
   read(datamap_unit, *) rangex
   read(datamap_unit, *) rangey
   read(datamap_unit, *) zero_posx
   read(datamap_unit, *) zero_posy
   read(datamap_unit, *) cfl
   close(datamap_unit)

   if (cpu_id .eq. 0) then
      write(*, *)
      write(*, *) 'Data space:  ', trim(directory) // trim(datavis_path)
      write(*, *)
      write(*, "(a15,i5)") ' x grid points=', (nx_total - ng * 2)
      write(*, "(a15,i5)") ' y grid points=', (ny_total - ng * 2)
      write(*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  readdata
!----------------------------------------------------------------------------------------------------------------------------------|
   open(datamap_unit, file = trim(datamap_path) // 'nt.mapout', form = 'unformatted')
   read(datamap_unit) nt
   close(datamap_unit)

   nx_lv = (nx_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng
   ny_lv = (ny_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng

   allocate(t(nt), dt(nt), p_pos(n_dims), n_blocks_all(nt), level_len(n_levels),                                                   &
      x(nx_block), y(ny_block), xh(nx_block), yh(ny_block), x_lv(nx_total), y_lv(ny_total),                                        &
      ro(nx_block, ny_block), mx(nx_block, ny_block), my(nx_block, ny_block), mz(nx_block, ny_block),                              &
      bx(nx_block, ny_block), by(nx_block, ny_block), bz(nx_block, ny_block), en(nx_block, ny_block),                              &
      gp(nx_block, ny_block), az(nx_block, ny_block), gx(nx_block, ny_block), gy(nx_block, ny_block),                              &
      roh(nx_block, ny_block), gph(nx_block, ny_block), mph(nx_block, ny_block),                                                   &
      teh(nx_block, ny_block), betah(nx_block, ny_block), azh(nx_block, ny_block),                                                 &
      vxh(nx_block, ny_block), vyh(nx_block, ny_block), vzh(nx_block, ny_block),                                                   &
      bxh(nx_block, ny_block), byh(nx_block, ny_block), bzh(nx_block, ny_block),                                                   &
      jxh(nx_block, ny_block), jyh(nx_block, ny_block), jzh(nx_block, ny_block),                                                   &
      fgxh(nx_block, ny_block), fgyh(nx_block, ny_block),                                                                          &
      fggxh(nx_block, ny_block), fggyh(nx_block, ny_block),                                                                        &
      fmgxh(nx_block, ny_block), fmgyh(nx_block, ny_block),                                                                        &
      fmtxh(nx_block, ny_block), fmtyh(nx_block, ny_block), fmtzh(nx_block, ny_block),                                             &
      ro_lv(nx_lv, ny_lv), gp_lv(nx_lv, ny_lv), mp_lv(nx_lv, ny_lv),                                                               &
      te_lv(nx_lv, ny_lv), beta_lv(nx_lv, ny_lv), az_lv(nx_lv, ny_lv),                                                             &
      vx_lv(nx_lv, ny_lv), vy_lv(nx_lv, ny_lv), vz_lv(nx_lv, ny_lv),                                                               &
      bx_lv(nx_lv, ny_lv), by_lv(nx_lv, ny_lv), bz_lv(nx_lv, ny_lv),                                                               &
      jx_lv(nx_lv, ny_lv), jy_lv(nx_lv, ny_lv), jz_lv(nx_lv, ny_lv),                                                               &
      fgx_lv(nx_lv, ny_lv), fgy_lv(nx_lv, ny_lv),                                                                                  &
      fggx_lv(nx_lv, ny_lv), fggy_lv(nx_lv, ny_lv),                                                                                &
      fmgx_lv(nx_lv, ny_lv), fmgy_lv(nx_lv, ny_lv),                                                                                &
      fmtx_lv(nx_lv, ny_lv), fmty_lv(nx_lv, ny_lv), fmtz_lv(nx_lv, ny_lv))

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   var_min = 1.0e10
   var_max = -1.0e10
   var_min_g = 1.0e10
   var_max_g = -1.0e10

   output_time_start = max(0, output_time_start) + 1
   output_time_end = min(nt - 1, output_time_end) + 1

   n_out = 0

   call mpi_barrier(mpi_comm_world, error)

   do ti = 1, nt
      write(ti_char, '(i4.4)') ti
      open(datamap_unit, file = trim(datamap_path) // 'info.t' // ti_char // '.mapout', form = 'unformatted')
      read(datamap_unit) t(ti)
      read(datamap_unit) dt(ti)
      close(datamap_unit)
   enddo

   do ti = output_time_start, output_time_end

      if (ti .ne. output_time_start + cpu_id + cpu_n * n_out) cycle
      n_out = n_out + 1

      write(ti_char, '(i4.4)') ti

!----------------------------------------------------------------------------------------------------------------------------------|
!     read t, dt
!----------------------------------------------------------------------------------------------------------------------------------|
!      open(datamap_unit, file = trim(datamap_path) // 'info.t' // ti_char // '.mapout', form = 'unformatted')
!      read(datamap_unit) t(ti)
!      read(datamap_unit) dt(ti)
!      close(datamap_unit)

      datavis_name = trim (datavis_path) // 'data' // ti_char
      write(*, *) ' Time = ', t(ti)
      i_block = 0
      nc_lv = 0

      if (output_method .eq. 2) open(datavis_unit, file = trim(datavis_name) // '.mapvis', form = 'unformatted')

      if (output_method .eq. 3) then
         data_offset = 0
         if (format .eq. 0) then
            open(datavtu_unit, file = trim(datavis_name) // '.vtu', form = 'binary')
            write(datavtu_unit) '<?xml version="1.0"?>' // char(10)
            write(datavtu_unit) '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">' // char(10)
            write(datavtu_unit) '<UnstructuredGrid>' // char(10)
         else
            open(datavtu_unit, file = trim(datavis_name) // '.vtu')
            write(datavtu_unit, '(a)') '<?xml version="1.0"?>'
            write(datavtu_unit, '(a)') '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
            write(datavtu_unit, '(a)') '<UnstructuredGrid>'
         endif

         if (format .eq. 0) then
            do cpui = 0, cpu_number - 1
               write(cpui_char, '(i4.4)') cpui

               open(datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.cpu' // cpui_char // '.mapout',            &
                  form = 'unformatted')
               read(datamap_unit) n_blocks_all(ti)
               read(datamap_unit) level_len
               read(datamap_unit) n_blocks
               do ci = 1, n_blocks
                  read(datamap_unit) id
                  read(datamap_unit) lv, nc, p_pos
                  read(datamap_unit) dx
                  read(datamap_unit) dy

                  read(datamap_unit) x
                  read(datamap_unit) y
                  read(datamap_unit) ro
                  read(datamap_unit) mx
                  read(datamap_unit) my
                  read(datamap_unit) mz
                  read(datamap_unit) bx
                  read(datamap_unit) by
                  read(datamap_unit) bz
                  read(datamap_unit) en
                  read(datamap_unit) az
                  read(datamap_unit) gx
                  read(datamap_unit) gy
                  if (nc .ne. 0) cycle
                  call write_vtu_header(datavtu_unit, var_name,                                                                    &
                     nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, n_var, var_dim, format, data_offset)
               enddo
               close(datamap_unit)
            enddo
         endif
      endif

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</UnstructuredGrid>' // char(10)
            write(datavtu_unit) '<AppendedData encoding="raw">' // char(10)
            write(datavtu_unit) '_'
         endif
      endif

      do cpui = 0, cpu_number - 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        read n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
         write(cpui_char, '(i4.4)') cpui

         open(datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.cpu' // cpui_char // '.mapout',                  &
            form = 'unformatted')
         read(datamap_unit) n_blocks_all(ti)
         read(datamap_unit) level_len
         read(datamap_unit) n_blocks
         write(datavis_unit) n_blocks

         if (cpui .eq. 0) then
            if (output_method .eq. 1) then
               allocate(var(sum(var_dim), nx_block, ny_block, n_blocks_all(ti)),                                                   &
                        x_arr(nx_block, n_blocks_all(ti)), y_arr(ny_block, n_blocks_all(ti)))
            endif
            if (output_method .eq. 3) then
               allocate(var(sum(var_dim), nx_block, ny_block, 1),                                                                  &
                        x_arr(nx_block, 1), y_arr(ny_block, 1))
            endif
         endif

         do ci = 1, n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
!           readx, ro, mx, my, mz, bx, by, bz, en
!----------------------------------------------------------------------------------------------------------------------------------|
            read(datamap_unit) id
            read(datamap_unit) lv, nc, p_pos
            read(datamap_unit) dx
            read(datamap_unit) dy

            read(datamap_unit) x
            read(datamap_unit) y
            read(datamap_unit) ro
            read(datamap_unit) mx
            read(datamap_unit) my
            read(datamap_unit) mz
            read(datamap_unit) bx
            read(datamap_unit) by
            read(datamap_unit) bz
            read(datamap_unit) en
            read(datamap_unit) az
            read(datamap_unit) gx
            read(datamap_unit) gy

            i_e = nx_block - 1
            j_e = ny_block - 1

            gp = (en - (mx ** 2 + my ** 2 + mz ** 2) / ro / 2.0d0 - (bx ** 2 + by ** 2 + bz ** 2) / 2.0d0) * gmm1

            do i = 1, i_e
               xh(i) = (x(i) + x(i + 1)) / 2.0d0
            enddo

            do j = 1, j_e
               yh(j) = (y(j) + y(j + 1)) / 2.0d0
            enddo

            do j = 1, j_e
               jp1 = j + 1
               do i = 1, i_e
                  ip1 = i + 1
                  roh(i, j) = (ro(i, j) + ro(ip1, j) + ro(i, jp1) + ro(ip1, jp1)) / 4.0
                  gph(i, j) = (gp(i, j) + gp(ip1, j) + gp(i, jp1) + gp(ip1, jp1)) / 4.0
                  bxh(i, j) = (bx(i, j) + bx(ip1, j) + bx(i, jp1) + bx(ip1, jp1)) / 4.0
                  byh(i, j) = (by(i, j) + by(ip1, j) + by(i, jp1) + by(ip1, jp1)) / 4.0
                  bzh(i, j) = (bz(i, j) + bz(ip1, j) + bz(i, jp1) + bz(ip1, jp1)) / 4.0
                  vxh(i, j) = (mx(i, j) + mx(ip1, j) + mx(i, jp1) + mx(ip1, jp1)) / 4.0 / roh(i, j)
                  vyh(i, j) = (my(i, j) + my(ip1, j) + my(i, jp1) + my(ip1, jp1)) / 4.0 / roh(i, j)
                  vzh(i, j) = (mz(i, j) + mz(ip1, j) + mz(i, jp1) + mz(ip1, jp1)) / 4.0 / roh(i, j)
                  mph(i, j) = (bxh(i, j) * bxh(i, j) + byh(i, j) * byh(i, j) + bzh(i, j) * bzh(i, j)) / 2.0
                  teh(i, j) = gph(i, j) / roh(i, j)
                  betah(i, j) = gph(i, j) / (mph(i, j) + eps)
                  azh(i, j) = (az(i, j) + az(ip1, j) + az(i, jp1) + az(ip1, jp1)) / 4.0
                  fgxh(i, j) = (gx(i, j) + gx(ip1, j) + gx(i, jp1) + gx(ip1, jp1)) / 4.0 * roh(i, j)
                  fgyh(i, j) = (gy(i, j) + gy(ip1, j) + gy(i, jp1) + gy(ip1, jp1)) / 4.0 * roh(i, j)
               enddo
            enddo

            do j = 2, j_e - 1
               do i = 2, i_e - 1
                  jxh(i, j) = (bzh(i, j + 1) - bzh(i, j - 1)) / dy / 2.0
                  jyh(i, j) =-(bzh(i + 1, j) - bzh(i - 1, j)) / dx / 2.0
                  jzh(i, j) = (byh(i + 1, j) - byh(i - 1, j)) / dx / 2.0 -                                                         &
                              (bxh(i, j + 1) - bxh(i, j - 1)) / dy / 2.0
                  fggxh(i, j) = (gph(i + 1, j) - gph(i - 1, j)) / dx / 2.0
                  fggyh(i, j) = (gph(i, j + 1) - gph(i, j - 1)) / dy / 2.0
                  fmgxh(i, j) = (mph(i + 1, j) - mph(i - 1, j)) / dx / 2.0
                  fmgyh(i, j) = (mph(i, j + 1) - mph(i, j - 1)) / dy / 2.0
                  fmtxh(i, j) = bxh(i, j) * (bxh(i + 1, j) - bxh(i - 1, j)) / dx / 2.0 +                                           &
                                byh(i, j) * (bxh(i, j + 1) - bxh(i, j - 1)) / dy / 2.0
                  fmtyh(i, j) = bxh(i, j) * (byh(i + 1, j) - byh(i - 1, j)) / dx / 2.0 +                                           &
                                byh(i, j) * (byh(i, j + 1) - byh(i, j - 1)) / dy / 2.0
                  fmtzh(i, j) = bxh(i, j) * (bzh(i + 1, j) - bzh(i - 1, j)) / dx / 2.0 +                                           &
                                byh(i, j) * (bzh(i, j + 1) - bzh(i, j - 1)) / dy / 2.0
               enddo
            enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!           for idl base level
!----------------------------------------------------------------------------------------------------------------------------------|
            if (lv .eq. output_level) then
               dx_lv = dx
               dy_lv = dy
               nc_lv = nc_lv + 1

               nx_tmp = (p_pos(1) - 1) * nx_no_bnd + ng
               ny_tmp = (p_pos(2) - 1) * ny_no_bnd + ng

               do i = 0, nx_no_bnd
                  i_c = i + ng
                  i_p = i + nx_tmp
                  x_lv(i_p) = xh(i_c)
               enddo

               do j = 0, ny_no_bnd
                  j_c = j + ng
                  j_p = j + ny_tmp
                  y_lv(j_p) = yh(j_c)
               enddo

               do j = 0, ny_no_bnd
                  j_c = j + ng
                  j_p = j + ny_tmp
                  do i = 0, nx_no_bnd
                     i_c = i + ng
                     i_p = i + nx_tmp
                     ro_lv(i_p, j_p) = roh(i_c, j_c)
                     gp_lv(i_p, j_p) = gph(i_c, j_c)
                     mp_lv(i_p, j_p) = mph(i_c, j_c)
                     te_lv(i_p, j_p) = teh(i_c, j_c)
                     beta_lv(i_p, j_p) = betah(i_c, j_c)
                     az_lv(i_p, j_p) = azh(i_c, j_c)
                     vx_lv(i_p, j_p) = vxh(i_c, j_c)
                     vy_lv(i_p, j_p) = vyh(i_c, j_c)
                     vz_lv(i_p, j_p) = vzh(i_c, j_c)
                     bx_lv(i_p, j_p) = bxh(i_c, j_c)
                     by_lv(i_p, j_p) = byh(i_c, j_c)
                     bz_lv(i_p, j_p) = bzh(i_c, j_c)
                     jx_lv(i_p, j_p) = jxh(i_c, j_c)
                     jy_lv(i_p, j_p) = jyh(i_c, j_c)
                     jz_lv(i_p, j_p) = jzh(i_c, j_c)
                     fgx_lv(i_p, j_p) = fgxh(i_c, j_c)
                     fgy_lv(i_p, j_p) = fgyh(i_c, j_c)
                     fggx_lv(i_p, j_p) = fggxh(i_c, j_c)
                     fggy_lv(i_p, j_p) = fggyh(i_c, j_c)
                     fmgx_lv(i_p, j_p) = fmgxh(i_c, j_c)
                     fmgy_lv(i_p, j_p) = fmgyh(i_c, j_c)
                     fmtx_lv(i_p, j_p) = fmtxh(i_c, j_c)
                     fmty_lv(i_p, j_p) = fmtyh(i_c, j_c)
                     fmtz_lv(i_p, j_p) = fmtzh(i_c, j_c)
                  enddo
               enddo
            endif

            if (output_method .eq. 2) then
               call write_mapvis(lv, nc, dx, dy, xh, yh,                                                                           &
                  roh, gph, mph, teh, betah, azh, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                                     &
                  fgxh, fgyh, fggxh, fggyh, fmgxh, fmgyh, fmtxh, fmtyh, fmtzh,                                                     &
                  var_flag, nx_block, ny_block, ng, n_var_total, datavis_unit)
            endif

            if (nc .ne. 0) cycle

!----------------------------------------------------------------------------------------------------------------------------------|
!           calculate the maximum and minimum value
!----------------------------------------------------------------------------------------------------------------------------------|
            call extremum_calculate(var_min, var_max, roh, gph, mph, teh, betah, azh, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh, &
               fgxh, fgyh, fggxh, fggyh, fmgxh, fmgyh, fmtxh, fmtyh, fmtzh, var_flag,                                              &
               n_var, nx_block, ny_block, ng, n_var_total)

!----------------------------------------------------------------------------------------------------------------------------------|
!           var array and coordinate array
!----------------------------------------------------------------------------------------------------------------------------------|
            if (output_method .eq. 1) then
               i_block = i_block + 1
               call data_collect(var, roh, gph, mph, teh, betah, azh, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                 &
                  fgxh, fgyh, fggxh, fggyh, fmgxh, fmgyh, fmtxh, fmtyh, fmtzh, var_flag,                                           &
                  x_arr, y_arr, xh, yh, sum(var_dim), nx_block, ny_block, n_blocks_all(ti), n_var_total, i_block)
            endif
            if (output_method .eq. 3) then
               i_block = 1
               call data_collect(var, roh, gph, mph, teh, betah, azh, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                 &
                  fgxh, fgyh, fggxh, fggyh, fmgxh, fmgyh, fmtxh, fmtyh, fmtzh, var_flag,                                           &
                  x_arr, y_arr, xh, yh, sum(var_dim), nx_block, ny_block, i_block, n_var_total, i_block)
            endif

            if (output_method .eq. 3) then
               call write_vtu_data(datavtu_unit, var_name,                                                                         &
                  real(var(:, ng : nx_block - ng, ng : ny_block - ng, i_block)),                                                   &
                  real(x_arr(ng : nx_block - ng, i_block)), real(y_arr(ng : ny_block - ng, i_block)),                              &
                  nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, n_var, var_dim, format)
            endif

         enddo
         close(datamap_unit)
      enddo

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</AppendedData>' // char(10)
            write(datavtu_unit) '</VTKFile>'
            close(datavtu_unit)
         else
            write(datavtu_unit, '(a)') '</UnstructuredGrid>'
            write(datavtu_unit, '(a)') '</VTKFile>'
            close(datavtu_unit)
         endif
      endif
      if (output_method .eq. 2) then
         call write_mapvis(output_level, nc_lv, dx_lv, dy_lv, x_lv, y_lv,                                                          &
            ro_lv, gp_lv, mp_lv, te_lv, beta_lv, az_lv, vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv, jx_lv, jy_lv, jz_lv,             &
            fgx_lv, fgy_lv, fggx_lv, fggy_lv, fmgx_lv, fmgy_lv, fmtx_lv, fmty_lv, fmtz_lv,                                         &
            var_flag, nx_lv, ny_lv, ng, n_var_total, datavis_unit)
         close(datavis_unit)
      endif
      if (output_method .eq. 1) then
         call write_vtk_unstructured(trim(datavis_name) // '.vtk', var_name,                                                       &
            real(var(:, ng : nx_block - ng, ng : ny_block - ng, :)),                                                               &
            real(x_arr(ng : nx_block - ng, :)), real(y_arr(ng : ny_block - ng, :)),                                                &
            nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, n_blocks_all(ti), n_var, var_dim, format, header)
      endif

      if (output_method .eq. 3) deallocate(var, x_arr, y_arr)
      if (output_method .eq. 1) deallocate(var, x_arr, y_arr)
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!   write parameters
!----------------------------------------------------------------------------------------------------------------------------------|
   do i_var = 1, n_var
      call mpi_allreduce(var_min(i_var), var_min_g(i_var), 1, mpi_real4, mpi_min, mpi_comm_world, error)
      var_min(i_var) = var_min_g(i_var)
      call mpi_allreduce(var_max(i_var), var_max_g(i_var), 1, mpi_real4, mpi_max, mpi_comm_world, error)
      var_max(i_var) = var_max_g(i_var)
   enddo

   if (cpu_id .eq. 0) then
      open(datavis_unit, file = trim(datavis_path) // 'par.mapvis')
      write(datavis_unit, '(a)') trim(directory) // trim(datavis_path)
      write(datavis_unit, '(i5, t32, a)') cpu_number, 'cpu_number'
      write(datavis_unit, '(i5, t32, a)') n_dims, 'n_dims'
      write(datavis_unit, '(i5, t32, a)') sys, 'sys'
      write(datavis_unit, '(i5, t32, a)') ng, 'ng'
      write(datavis_unit, '(i5, t32, a)') ini_block_x, 'ini_block_x'
      write(datavis_unit, '(i5, t32, a)') ini_block_y, 'ini_block_y'
      write(datavis_unit, '(i5, t32, a)') nx_no_bnd, 'nx_no_bnd'
      write(datavis_unit, '(i5, t32, a)') nx_block, 'nx_block'
      write(datavis_unit, '(i5, t32, a)') nx_total, 'nx_total'
      write(datavis_unit, '(i5, t32, a)') ny_no_bnd, 'ny_no_bnd'
      write(datavis_unit, '(i5, t32, a)') ny_block, 'ny_block'
      write(datavis_unit, '(i5, t32, a)') ny_total, 'ny_total'
      write(datavis_unit, '(i5, t32, a)') ratio, 'ratio'
      write(datavis_unit, '(i5, t32, a)') n_levels, 'n_levels'
      write(datavis_unit, '(i5, t32, a)') div_clean_flag, 'div_clean_flag'
      write(datavis_unit, '(i5, t32, a)') output_real8_flag, 'output_real8_flag'
      write(datavis_unit, '(i5, t32, a)') output_level, 'output_level'
      write(datavis_unit, '(es15.6, t32, a)') gm, 'gm'
      write(datavis_unit, '(es15.6, t32, a)') rangex, 'rangex'
      write(datavis_unit, '(es15.6, t32, a)') rangey, 'rangey'
      write(datavis_unit, '(es15.6, t32, a)') zero_posx, 'zero_posx'
      write(datavis_unit, '(es15.6, t32, a)') zero_posy, 'zero_posy'
      write(datavis_unit, '(es15.6, t32, a)') cfl, 'cfl'
      write(datavis_unit, '(i5, t32, a)') nt, 'nt'
      do ti = 1, nt
         write(datavis_unit, '(es15.6, t32, a, i3, a)') t(ti), 't (ti = ', ti, ')'
      enddo
      do ti = 1, nt
         write(datavis_unit, '(es15.6, t32, a, i3, a)') dt(ti), 'dt (ti =', ti, ')'
      enddo
      write(datavis_unit, '(i5, t32, a)') n_var, 'n_var'
      write(datavis_unit, '(t4, a5, t32, a, i3, a)') (var_name(i_var), 'var_name (i_var = ', i_var, ')', i_var = 1, n_var)
      write(datavis_unit, '(es15.6, t32, a, a5, a)') (var_max(i_var), 'var_max (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
      write(datavis_unit, '(es15.6, t32, a, a5, a)') (var_min(i_var), 'var_min (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
      write(datavis_unit, '(i5, t32, a, a5, a)') (var_dim(i_var), 'var_dim (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
      close(datavis_unit)
   endif

   deallocate(t, dt, p_pos, n_blocks_all, level_len, x, y, xh, yh, x_lv, y_lv,                                                     &
      ro, mx, my, mz, bx, by, bz, en, gp, az, gx, gy, roh, gph, mph, teh, betah, azh, vxh, vyh, vzh,                               &
      bxh, byh, bzh, jxh, jyh, jzh, fgxh, fgyh, fggxh, fggyh, fmgxh, fmgyh, fmtxh, fmtyh, fmtzh,                                   &
      ro_lv, gp_lv, mp_lv, te_lv, beta_lv, az_lv,  vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv, jx_lv, jy_lv, jz_lv,                  &
      fgx_lv, fgy_lv, fggx_lv, fggy_lv, fmgx_lv, fmgy_lv, fmtx_lv, fmty_lv, fmtz_lv, var_min, var_max, var_min_g, var_max_g)

   call mpi_barrier(mpi_comm_world, error)

   if (cpu_id .eq. 0) then
      write(*,*)
      write(*,*) ' readdata complete!'
      write(*,*)
   endif

   call mpi_finalize(error)

!----------------------------------------------------------------------------------------------------------------------------------|
   stop
end program data_manage
